Finite-size scaling for non-linear rheology of fluids confined in a small space 
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We perform molecular dynamics simulations in order to examine the rheological transition of 
fluids confined in a small space. By performing finite-size scaling analysis, we demonstrate that this 
rheological transition results from the competition between the system size and the length scale of 
cooperative particle motion. 
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When fluids are confined in a small space, their rheo- 
logical property changes drastically from that of the flu- 
ids in the bulk. A typical example of this rheological 
transition is observed in a fluid confined between solid 
walls. When the distance between the walls is so large 
that the fluid is considered in the bulk, it behaves as 
a normal Newtonian fluid. On the other hand, when 
the distance becomes equivalent to that of approximately 
eight molecular layers, the fluid exhibits an increase in its 
viscosity, the shear-thinning behavior and the appear- 
ance of the yield stress as well as the enhancement of the 
relaxation time [l| . This rheological transition has been 
extensively studied by experiments 0, H & n d simulations 

assaa. 

Some researchers have conjectured that this rheologi- 
cal transition results from structural transition 0, @, 0] • 
However, in some simulations of fluids confined between 
walls, the rheological transition is observed without any 
clear structural transition [4|. Hence, we consider that 
there is another mechanism responsible for this rheolog- 
ical transition. 

Let us recall that glassy materials, such as dense col- 
loidal suspensions and super-cooled liquids, also display 
rheological transition similar to that observed in fluids 
confined in a small space when the temperature decreases 
or the density increases 0, EH EH- From this fact, we 
conjecture that fluids confined in a small space and glassy 
materials have common features [H, 0] • 

In particular, in glassy materials, it has been observed 
that the length scale of cooperative particle motion, 
which is called the dynamical correlation length, becomes 
comparable with the system size near glass transitions 
[12, EH- Therefore, in the case of fluids confined in a 
small space, the dynamical correlation length might be 
comparable with the system size when the size is suffi- 
ciently small. Our conjecture suggests that the compe- 
tition between the dynamical correlation length and the 
system size is an origin of the rheological transition. 

In order to test whether the rheological transition is 
actually caused by the conjectured effect, we investigate 
a model that does not exhibit any structural transition. 
By performing numerical simulations with a finite-size 
scaling analysis [14| , we demonstrate that the rheological 
transition occurs because of the competition between the 



system size and the dynamical correlation length. 

Model: The system is a 80:20 mixture of Lennard- 
Jones particles of types A and B in a two-dimensional 
square box, with an interaction potential 



V a p{r) = 4e Q/3 



(T £ ) 12 -( f r)' 



in 



where a and j3 refer to the two species A and B. In order 
to restrict crystallization, we investigate the binary mix- 
ture. The particles have equal masses, and the interac- 
tion parameters are 6ab = 1-5£ab, £bb = 0.5£aa, cab = 
0.88oaa, and ctbb = 0.8ctaa- The length, energy, and 
time units are expressed in the standard Lennard- Jones 
units: oaa (particle diameter), €aa (interaction energy), 
and To = (mACT^/eAA) 1 ^ 2 , where toa is the particle 
mass. The subscript A refers to the major species. We 
fix the total number density p = N/L 2 = 1.13, where L 
is the system size and N is the number of the particles. 

The evolution equations of the z-th particle's position 
Ti = (xi,yi) and velocity Vi = (v,i,Vi) are numerically 
integrated by using a leapfrog algorithm with a time step 
5t = 0.005. Here, in order to maintain the temperature 
T, we use the Nose- Hoover algorithm [IB]] . In order to 
realize the shear flow (jy, 0) without structural ordering 
near the boundaries, Lees-Edwards boundary conditions 
are imposed [HI]. Here, 7 is the shear rate. 

Equilibrium Properties: We first calculate the pair 
distribution function defined as 



9( 



1,L ■ 



(2) 



Here, (-) T L represents the ensemble average under the 
condition that temperature T, the shear rate 7, and the 
system size L are provided. In Fig. [H we show the 
equilibrium pair distribution functions for several values 
of L when T = 0.8 and 7 = 0. It can be seen that the 
structural transition does not occur when the system size 
is decreased. 

Next, we study the dynamical correlation length £(T), 
following Ref . [H, EH . For this purpose, we define the 
local indicator of configuration change as 



F(r, t) = J2 5 ( r - ri(0))W(n(t) - r,(0)), 



(3) 
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FIG. 1: Pair distribution function g(r) for several values of 
L when T = 0.8 and 7 = 0. 



FIG. 2: x(ti T) as a function of t for several values of T when 
L = 37.6. 



with 



W(r) 



1 (\r\ < a), 
(otherwise). 



(4) 



Here, we choose the value of a as 0.3, which is slightly 
larger than the square root of the plateau value of the 
mean square displacement, as discussed in Ref. fl2| . 
Note that (F(r,t)) is akin to the intermediate scattering 
function when the system is homogeneous. The correla- 
tion function C(r,t,T) of F(r,t) is defined as 



V 2 



C(r, t,T) = — (F(r, t)F(0, t)) T ^ L ~ (F(r, t))^ L . 

(5) 

From C(r,t,T), the amplitude of the correlation is de- 
fined as 



x(t,T) 



k 



l — J drC(r,t,T). 



(6) 



In Fig. [51 we show the time dependence of T) for 
several values of T when L — 37.6 and 7 = 0. x(t> T) 
has a maximum value at a time i max . As the tempera- 
ture decreases, i max and the maximum value x(^max,T) 
increases. 

In order to extract £(T), we define the Fourier trans- 
formation S(k,t,T) of C(r,t,T) as 

S(k,t,T) = J drC(r,t,T)exp{-ik-r), (7) 
S(0,t,T) = k B T x (t,T). (8) 



where 



From S(t,k,T), we define S max (k,T) = S(k,t max ,T). By 
assuming the functional form of S max (k,T) as 



Smax{k, T) — 



ffmax(0,r) 



(9) 



we estimate £(T) for the case that 7 = and L = 37.6. In 
Fig. S max (k, T)/5 max (0, T) is displayed as a function 
of fc£(T). In Fig. E2 we show £(T) as a function of T. 
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FIG. 3: S max (k,T)/S max (0,T) as a function of fc£(T) for 
several values of T when L — 37.6. 



Rheological Properties: As a typical rheological prop- 
erty, we investigate the viscosity. In our system, the shear 
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FIG. 4: f (T) as a function of T for L = 37.6. 
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stress a X y(T,^, L) is expressed as 



a xy (T,j,L) = — ( 



rriiUiVi 



E 



XjjVij dV(r i3 ) 



1 T,7,L 

(10) 

where x t j = xi - Xj, y t j — y t - yj, and r io = \n - Tj\ 
[lit . From a xy (T,j,L) and 7, the viscosity 7/(T, 7,L) is 
defined as 



a xy {T,j,L) 

1 



(11) 
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In Fig. we show n(T, 7, L) as a function of 7 for several FIG. 6: Relation between the viscosity y(T, 7, L) and the 
values of T when L = 37.6. For T = 1.0, the viscosity shear rate 7 for several values of L when T = 0.6 and 0.9. 
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FIG. 5: Relation between y and 7 for several values of T 
when L = 37.6. 

is independent of the shear rate when 7 < 0.001. This 
regime is called a Newtonian regime. As the shear rate 
increases, the viscosity decreases as 77 ~ 7~ Q - When the 
temperature decreases, the Newtonian regime becomes 
narrower and finally disappears. The exponent a = 2/3 
is observed in experiments and simulations [3] of thin 
liquid films. Furthermore, this exponent is observed in 
simulations of glassy materials [ill ]. 

Now, we investigate the system size dependence of the 
rheological property. As examples, in Fig. [6l we show 
rj{T, 7, L) as a function of 7 for several values of L when 
T = 0.6 and 0.9. At T = 0.6, 77 changes slightly. On 
the other hand, when T = 0.9, 77 strongly depends on L 
in the low shear rate regime as in experiments 0, Q and 
simulations 0, H, H, 0, H| of a thin liquid film. 

Finite-size scaling: We perform a finite-size scaling 
analysis in a similar manner to that used for critical phe- 
nomena [3]. We first assume the scaling form of the 
viscosity as 

T7(T, 7 ,L) =7 /0 (T,i)G( 7 -%(T 1 ,L)i) , (12) 

where 770 (T,L) is the Newtonian viscosity defined as 
r]o(T,L) = lim 7 _>o 7y(T, 7, L). G(x) is a function having 



the properties 

MmG(x) = 1, (13) 

X— »0 

lim G(x) - x-i. (14) 

(15) 

It must be noted that this scaling form can be derived 
from a theory for the rheological property in glassy mate- 
rials [TH. Here, we only consider the case where 77o(7\ L) 
has a finite value. Next, we assume the functional form 
of ?/o (T, L) by £(T) and L as 

m (T,L) = ru(T)H(t(T)/L), (16) 

where r\i (T) is the viscosity in the thermodynamical limit 
L — > 00 and H(x) is the function having the property 

WmH(x) = 1. (17) 

Then, we test the scaling assumptions described by 
Eqs. (H2J) and ([16]). We show rj(T, 7, L)/rj (T, L) as a 
function of 7 ■ t]q(T, L) 3 / 2 in Fig. [7] In addition, we show 
77o(T, L)/ru(T) as a function of £{T)/L in Fig. E These 
figures clearly indicate that the scaling assumptions de- 
scribed by Eqs. (fT2]) and (fT6|) are plausible. In addition, 
the viscosity 77/ (T) in the thermodynamical limit is shown 
in the inset of Figs. [U From these results, we conclude 
that the size dependence of the rheological property in 
our system results from the competition between L and 

an 

Conclusion and Discussion: We performed the 
molecular dynamics simulations of the fluid confined in 
a small space with the Lees-Edwards boundary condi- 
tions. In these simulations, we observed the strong size 
dependence of the rheological property. In addition, we 
numerically calculated the correlation length £ (T) of the 
local configuration change. From the numerical confirma- 
tion of the scaling relations given by Eqs. (fl"2"|) and (|16p . 
it has been concluded that the strong size dependence 
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local indicator of the local configuration change, which is 
the dynamical quantity, is used in our finite-size scaling 
analysis. This choice is natural because the viscosity is 
approximately calculated from the time correlation func- 
tion of the local configuration change [17] ■ On the other 
hand, when we recall the fact that the viscosity is calcu- 
lated from the shear stress, which is the static quantity, it 
is also expected that the correlation of a static quantity 
becomes involved in the scaling analysis. However, as far 
as we have studied, we could not find such a quantity. 
This problem will be studied in future. 



FIG. 7: r)(T, 7, L) /rjo (T, L) as a function of 7 • 
various values of T and system size L. 
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FIG. 8: Main: 770(7, L)/rn{T) as a function of £{T)/L. Inset 
: rji (T) as a function of T. 



of the rheological property is caused by the competition 
between the system size and the correlation length. 

Although we have shown that the proposed mechanism 
affects the rheological property in our system without 
wall boundaries, it is not clear whether it would be re- 
ally effective in experiments of a thin fluid film with wall 
boundaries. In such a situation, other effects such as the 
structural transition might be important for determining 
the rheological properties. 

In Ref. [18|, it is shown that the density relaxation 
time of glassy materials increases when the system size 
becomes comparable with a dynamical correlation length. 
With the assumption that the relaxation time is propor- 
tional to the viscosity , the result in Ref. [l8| is similar 
to the result described by Eq. (1 16(1 and Fig. [U However, 
we wish to note that in Ref. [lq . the rheological prop- 
erty was not studied, and the finite size scaling was not 
performed. 

Finally, let us remind that the correlation length of the 



We thank S. Sasa and K. Hukushima for their valuable 
advice, and H. Tasaki, H. Hayakawa, and T. Hatano for 
their useful comments. 



[1] 

[2] 

[3 
[5 

[7; 

[8 
[9 

[10; 

[11 

[12. 

[13 
[14 

[15; 

[16 
[17 

[is; 

[19 



Electronic address: otsuki@jiro.c.u-tokyo. ac.jp 
Y-Z. Hu and S. Granick, Tribol. Lett. 5, 81 (1998). 
H-W. Hu, G. A. Carson, and S. Granick, Phys. Rev. Lett. 
66, 2758 (1991). 

S. Yamada, G. Nakamura, Y. Hanada, and T. Amiya, 
Tribol. Lett. 15, 83 (2003). 

P. A. Thompson, G. S. Grest, and M. O. Robbins, Phys. 
Rev. Lett. 68, 3448 (1992). 

A. Jabbarzadeh, P. Harrowell, and R. I. Tanner, Phys. 
Rev. Lett. 96, 206102 (2006). 

J. Klein and E. Kumacheva, J. Chem. Phys. 108, 6996 
(1998). 

S. T. Cui, P. T. Cummings, and H. D. Cochran, J. Chem. 
Phys. 114, 7189 (2001). 

S. T. Cui, C. McCabe, P. T. Cummings, and H. D. 
Cochran, J. Chem. Phys. 118, 8941 (2003). 
W. B. Russel, D. A. Saville, and W. R. Schowalter, Col- 
loidal Dispersions (Cambridge University Press, New 
York). 

R. G. Larson, The Structure and Rheology of Complex 
Fluids (Oxford University Press, New York, 1999). 
L. Berthier, and J.-L. Barrat, J. Chem. Phys. 116, 6228 
(2002). 

N. Lacevic, F. W. Starr, T. B. Schroder, and S. C. 

Glotzer, J. Chem. Phys. 119, 7372 (2003). 

L. Berthier, Phys. Rev. E 69, 020201(R) (2004). 

N. Goldenfeld, Lectures on Phase Transitions and the 

Renormalizatioin Group (Addison- Wesley, New York, 

1992). 

D. J. Evans and G. Morris, Statistical Mechanics of 
Nonequilibrium Liquids, (Academic, London, 1990). 
M. Otsuki and S. Sasa, J. Stat. Mech., L10004 (2006). 
K. Miyazaki, D. R. Reichman, and R. Yamamoto, Phys. 
Rev. E 70, 011501 (2004). 

K. Kim, and R. Yamamoto, Phys. Rev. E 61, R41 (2000). 
R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 
(1998). 



